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By determining the type of all stationary points of the Gibbs free energy functional for layered 
superconductors in parallel magnetic fields, we establish the classification of all solutions to coupled 
static sine-Gordon equations for the phase differences with respect to their stability. We prove that 
the only minimizers of the free energy are the Meissner solution (the "vacuum" state) and soliton 
vortex-plane soluti ons [S. V. Kuplevakhsky, Phys. Rev. B 60, 7496 (1999); ibid. 63, 054508 (2001); 
sond-mat /020229^ | . They are the actual equilibrium field configurations. We present a topolog- 
ical classification of these solutions. In contrast, previously proposed non-soliton configurations 
("isolated fluxons", "triangular Josephson-vortex lattices", etc.) are absolutely unstable and unob- 
servable: They are nothing but saddle points of the Gibbs free-energy functional and are not even 
stationary points of the Helmholtz free-energy functional (obtained from the former by a Legendre 
transformation). (Physically, non-soliton configurations violate conservation laws for the current 
and the flux.) The obtained results allow us to explain dynamic stability of vortex planes, noticed 
in numerical simulations, and to provide a unified interpretation of the available experimental data. 
We hope that the paper will stimulate interest in the subject of specialists in different fields of 
physics and in applied mathematics. 
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I. INTRODUCTION 



In this paper, we present the solution of the problem of equilibrium vortex structure in layered superconductors and 
stacked Josephson junctions in the presence of a parallel, static, homogeneous external magnetic field H and provide 
a unified interpretation of the available experimental data. Our approach consists in a rigorous mathematical analysis 
of the stability of all types of flux confisuxations, proposed in the literature, by means of exact variational methods 
for microscopic free-energy functionals Ecra 

Within the framework of these methods, we have previously obtained a complete classification of all possible static 
soliton solutions to coupled sine-Gordon (SG) equations.Jjpr phase differences both in infinite [N = oo, N is, the number 

ijiinite (N < oo) Josephson-junction stackaa for H > 0. 
i'B that topological solitons in nonlinear field theories are 



of superconducting(S) layers] layered superconductors 
Based on the fundamental argument of soliton physics 

minimizers of the energy functionals (frcc-cncrgy functionals in our case),IL2l we have identified these solutions with 
equilibrium Josephson-vortex configurations. Their magnetic field has symmetry typical of plane defects, hence the 
term "vortex planes". Physically, a vortex plane can be regarded as a bound state of interlayer vortices (one vortex 
per each insulating layer in the plane). In contrast to a deep-rooted belieftil in an "analogy" with Abrikosov vortices 
in continuum type-II superconductors, the SG equations for H > do not admit static soliton solutions that can be 
identified with an " isolated Josephson vortex" or a " triangular Josephson-vortex lattice" . 

Unfortunately, a wide-spread misunderstanding of the fact that equilibrium Josephson vortices are nothing but 
static soliton solutions to the SG equations incurred a misunderstanding of the exact mathematical results of Refs. 



1,2M In the critical comment, Ref. 12], V. M. Krasnov and L. N. Bulaevskii "disprove" the conclusions of Refs. 
1,2 by claiming (in contradiction to the fundamentals of soliton physics) that vortex-plane solitons "maximize the 
free enepay" £3 According to Ref. [l^, the "instability" of vortex planes is "similar to the instability of the laminar 
solutionis for type-II superconductors" . [Alternating superconducting and normal layers, envisaged by the laminar 
model,EJ have nothing to do with soliton physics and do not possess the property of topological stability: see a 
proof in subsection V.B of the present paper.] V. M. Krasnov and L. N. Bulaevskii insist on hypothetical "isolated 
fluxons", allegedly, having "lower energy" for H > and characterized by a "much smaller length scale". No exact 
mathematical definition of the "isolated fluxons" is given. We have not found any definition in the original papers by 
the critics of vortex planes, either: For example, it is claimed in Refs. |ll|[l6] that Josephson vortices "do not exist" in 
single Josephson junctions withJ^F 2Aj (W is the junction width, Aj is the Josephson length). However, an exact, 
closed-form analytical solutioroO to the single static SG equation, appropriate for this case, clearly demonstrates the 
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existence of phase-difference solitons for arbitrarily small W, provided the external field H is sufficiently high. [As 
shown in Refs. |l|,^, exactly these solitons account for the well-known Fraunhofer pattern of the critical Josephson 
current I c (H) for W < 2Xj.} 

Furthermore, the manuscript Ref. |[ although submitted twice to Physical Review B (November 2000, August 
200f), is still not published. In particular-aiLe of the referees disputed the conclusions of Ref. |[ because, in his 
opinion, the soliton boundary conditionsSO'LrH'EI employed therein " overdetermined" the problem of the classification 
of equilibrium Josephson-vortex configurations. He argued that certain numerical simulations for the SG equations 
had demonstrated, aside from vortex-plane solitons, the existence of " single- vortex" solutions. According to the 
referee, these solutions had "lower free energy" than the vortex planes for given H. As in Ref. [L2J, no exact definition 
of such solutions was given. 

It should be emphasized that the idea of an "analogy" between the Josephson-vortex structure in layered super- 
conductors and the Abrikospv.- vertex structure in continuum type-II superconductors was not supported in Ref. [If] 
and subsequent publication j 18 H 19 H 2C l i P 1 l by any serious mathematical arguments. EIBE3 Neither was it confirmed by direct, 
experimental observations of the equilibrium Josephson-vortex structure in artificial stacked junctions at H > OEj-EII 
Unfortunately, most theoretical efforts were constrained by the idea of an " apft lQSy i'i j . he nce the use of mathemati- 
cally ill-formulated methods, such as, e.g., a " continuum- limit approximation" .yllaSESEII For instance, the exact SG 
equations for the phase differences were not even derived in Refs. [Tl]Jl8|, concerned with a "single Josephson vortex" 
at H > 0. (As shown by Faridp3 equations of Ref. |f8| have no physical solution.) Since the problem of the stability 
of the proposed "vortex configurations" (i.e., whether they are actual points of minima of the free-energy functionals) 
required the use of rigorous mathematical methods, it was not ei££ii_posed in Refs. [if 18,lS|, pO[pT[ 

Concerning numerical simulations for the static SG equationsjijo there is an unjustified tendency to identify any 
kink- type feature of the phase difference with a " Josephson vortex" , without any analysis of its stability. In contrast 
to the exact analytical methods of Refs. (jJHIj, the numerical approach does not provide any means to-establish a 
full set of necessary and sufficient conditions of the minimum of the free-energy functionals. Typically,c3 numerical 
simulations start with an incorrectly formulated (both mathematically and physically) boundary value problem that 
does not meet the criterion of uniqueness. E£l 

To close the issue of the equilibrium Josephson-vortex structure in layered superconductors and stacked, iujaci. ions, 
we determine analytically (by means of exact methods of the calculus of variation* ?3 and soliton physicsoou'ElBEa) the 
type of all stationary points of the exact microscopic Gibbs free-energy functional,a generating the static SG equations. 
Our consideration applies to an arbitrary number of superconducting layers N, including the cases N = 2 (a single 
junction) and N — > oo (an infinite layered superconductor). As a result, we obtain a complete classification of all 
nontrivial solutions to the SG equations, considered in the literature (both analytically and numerically), with respect 
to their stability. As could be expected from the general arguments of soliton physics, the only minimizfirs. of the free- 
energy functional are the Meissner solution (the "vacuum" state) and soliton vortex-plane solutions .ana The latter 
represent the actual equilibrium Josephson-vortex configurations for H > 0. In contrast, non-soliton configurations 
(e.g., "isolated fluxons", "triangular Josephson-vortex lattices", etc.) are absolutely unstable: They are nothing 
but saddle points of the Gibbs free-energy functional and are not even stationary points of the Helmholtz free- 
energy functional (obtained from the former by a Legendre transformation) . [Physically, non-soliton solutions violate 
conservation laws for the current and the flux, as was first noticed in the case of an infinite layered superconductor 
(N = oo) in Refs. @,|] 

In section II, we clarify a relationship between the correct formulation of the boundary value problem to the SG 
equations and a full set of necessary and sufficient conditions of the minimum of the Gibbs free-energy functional. 
The proof of the stability of the Meissner solution and vortex-plane solitons is given in sections III, IV. In section III, 
we establish that the sufficient condition of the minimum of the Gibbs free-energy functional consists in the vanishing 
of the surface variation of the corresponding Helmholtz free-energy functional, which, in turn, yields conservation 
laws for the flux and the intralayer current. In subsection IV. A, we derive the soliton boundary conditions directly 
from the conservation law for the flux, which provides the sought proof of the stability of the Meissner solution and 
vortex-plane solitons. The main physical and mathematical properties of these solutions are discussed in subsection 
IV. B. I^puJyection IV. C, we analyze the obtained results from the point of view of general theory .of. topological 
defectseoQoEoEj and explain dynamic stability of vortex planes, established in numerical simulationsBiH'Ej 

Some important physical and mathematical issues, related to the main results of the paper, are discussed in 
section V. In subsection V.A, we present a rigorous analytical description of unstable solutions ("isolated fluxons", 
"triangular Josephson-vortex lattices" , etc.), proposed in previous theoretical publications and numerical simulations. 
In subsection V.B, we draw a comparison with the Abrikosov vortices in type-II superconductors. In subsection V.C, 
we analyze the available experimental data from the point of view of the stable vortex-plane configurations. 

In section VI, we summarize the obtained results, systematize our criticism of previous approaches and make some 
concluding remarks. In Appendix A, we give a list of mathematical formulas, relevant to the main text. In Appendix 
B, we establish a relationship to the variational principle of Refs. np for infinite layered superconductors. 
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Throughout the paper, we adhere to the dimensionless notation of Ref. |3|. The geometry of the problem is that 
of figures 1, 2 in Ref. [|: The superconductor occupies the region [0 < x < N — 1] x [— L < y < L] x (— oo < z < oo), 
where 2L — W; the layering axis (the c-axis) is x; the axis y is along the layers; the external magnetic field is along 
the axis z: H = (0,0, H > 0). The phase difference between two successive S-layers is denoted as <p n = ip n — ip n -i 
(</>o = 4>n = 0), with n = 0, 1, . . . , N — 1 being the S-layer number. The c-axis external current is not considered: 
1 = 0. 



II. THE FORMULATION OF THE PROBLEM 

In Ref. ^, we have derived the coupled, static SG equations for the phase differences (f>„ by minimizing the exact 
microscopic Gibbs free-energy functionalcl [f n , 4> n , A; H] with respect to /„ and A (/„ is the reduced modulus of 
the order parameter in the n-th S-layer, A is the vector potential). In the limit r (T) -C 1, H <C H c2 [r (T) is 
the parameter of the interlayer coupling, H c2 is the upper critical field], when f n = 1, the SG equations appear as 
solubility conditions for the Maxwell equations in the gauge 

A = [0,A(x,y),Q]. (I) 

The SG equations readJ§ 

^T^ = ^E G_1 ( n ' TO ) sin ^(y)' n=l,...,N-l; (2) 

y m-l 

where G" 1 (n, m) is a Jacobian matrix with elements G _1 (n,n) = 2 + e 2 (n = 1, . . . , N — 1), G _1 (n + l,n) = 
G" 1 (n,n + l) = -1 (n= l,...,JV-2), and G' 1 (n,m) = for \n - m\ > 1. 

The requirement that the local field H n (y) (n — l<x<n,n = l,...,N—l)be equal to the applied one at y = ±L 
has led to the conditions 

[The condition ^ (-L) > merely reflects the fact that the local field is parallel to the applied one H > 0. ] From 
the requirement that the local field be equal to the applied one at x = 0, x = N — 1, we have obtained 

<i>n (y) = <f>N-n (y) , n= 1,...,N- 1. (5) 

Equations (||) and boundary conditions (||) are satisfied by functions of the type 

<t>n (y) = -<f>n (-y) + 2ttZ„, (6) 

where the constants Z n can be arbitrarily chosen from the set 0, ±1, ±2, ... As is pointed out in Ref. [| the fixation 
of the constants Z n requires imposition of boundary conditions on <p n at y — ±L. Based on general arguments of 
soliton physics that soliton solutions are minimizers of corresponding energy functionals, we have imposed in Ref. || 
standard soliton boundary conditions on <p n . 

As it has turned out, the only possible solutions, compatible with the requirement (|J), are the Meissner solution 
and the soliton vortex-plane solutions, for which (^) is satisfied automatically and 

Z n = Z n+1 =N V , n = 1, . . . , N - 2; N v = 0, 1, 2, . . . (7) 

in (fsf) , with N v = representing the topologically trivial Meissner solution. Given thatS 



N-l 



H n (y) =H[G (n, 1) + G (n, N - 1)] + j ^ G (n, m) (8) 

m=l " 

the field H = H, corresponding to a concrete configuration {</>„} with N v = N v , is determined by 

f (-L) = 2H. (9) 
ay 
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Note that the first, phase- independent, term in p is a contribution of the field penetrating through the boundaries 
x = 0, x = N — 1, and the second term is a contribution of the field penetrating through the boundaries y = ±L. 
The matrix G (n, m)is the inverse of G~ l (n, to): for its properties, see Appendix A. 

In contrast to the above exact variational method, numerical simulationsE£lc3 for (||) start with the imposition of 
the boundary conditions 



^ (±L) = 2H, 
dy 



(10) 



without any regard to sufficient conditions of the minimum of the Gibbs free-energy functional. Such an approach is 



based on an erroneous belief that all solutions to 
belief arc the naive " enersu-areuments" of Ref. 12 



(|10| ) minimize the free-energy functional. (Manifestations of this 
appealing to " a difference in the length scales" , and calculations 
by means of combinatoricsE9{£3 of the "number of quasi-equilibrium fluxon modes" J However, conditions (|Io| ) do not 
specify any boundary value problem for (^) : By virtue of the symmetry relations (q) , the imposition of the boundary 
condition on ^j- at y = — L automatically ensures the fulfillment of the same boundary condition at y = L, whereas 
the constants Z n remain undetermined. Thus, the "boundary value problem" Jl0| ) does not satisfy the criterion of 
uniqueness, l£I which is a sign of the presence of unphysical (i.e., unobservable) solutions. 

The existence of redundant solutions to (||), (10) is already clear for physical reasons: This "boundary value 
problem" does not take any account of the necessity to ensure the continuity of the local field at the boundaries x = 0, 
x = N — 1. To understand at a rigorous mathematical level where the unphysical solutions come from, we have to 



consider all stationary points of the generating Gibbs free-energy functional, rewritten via <f> n and .Q In this way, 
we will derive a full set of necessary and sufficient conditions of the minimum directly from the variational_urinciple, 
obtain an independent proof of the fact that the Meissner solution and the soliton vortex-plane solutionsElau are the 
unique minimizers of the problem and establish the character of the instability of unphysical non-soli±on solutions. In 
our consideration, we will employ the first integral of (||) that, taking account of (0), has the formO 



N-l 



E 



N-l N-l 



<" s <t>n (y) + y ^ }_^G (n, TO 

n— 1 ri— 1 m— 1 



dc/) n (y) dcj) m (y) 
dy dy 



2H 2 
~W 



N-l 



(N-l)+ 5]cos0 n (£). 



(11) 



where H s is the superheating (penetration) field of a scmiinfinitc (0 < y < +oo) Josephson-junction stack, given by 
Eq. (|A8|). 



III. THE NECESSARY AND SUFFICIENT CONDITIONS OF THE MINIMUM OF THE GIBBS 

FREE-ENERGY FUNCTIONAL 



The generating Gibbs free-energy functional for the SG equations (|2j) has the for 

JV-l 



mi 



dy ' 



F 



dy 



H 



AHr (T) $1 



■» (L) - K {-L) 
2ir 



(12) 




N-i h. 
+ E / dy t 1 ~ cos ^« 

n=l J L 



(13) 
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n[l-G(n,l)-G(n,N -I)] 



NW 



Note that the functional ( |12[) is measured from the condensation energy fto = — 
external field in the absence of the sample is subtracted.^ 

Our task is to establish the necessary and sufficient conditions of the minimum of ( [l2|) 



(14) 

Moreover, the energy of the 



Aft 



dy 



= ft 



4>n + Hn, ~7~~ + 



dy dy 



- ft 


9ni —j—i H 




dy 



> 



(15) 



(16) 



First, we observe that, in contrast to (|12[), the functional (13) is positive definite 

dy 

since the matrix G (n, m) is positive definite: see Appendix A. The absolute minimum of (|l3| ) is achieved for H = 0, 
4> n = (n = 1, 2, . . . , N — 1). Hence, the functional necessarily has minima for any H > 0. The first variations 
of (113), (Eah are 



5ft 



dy 



= SF 



dy 



N-l 



- 4Hr (T) $ < 



n=l 



S(j> n (L) - 5<j> n (-L) 
2ir 



(17) 



SF 



dy 



H 



N-l u 

(t) E J d y 



N-l 



sin 4> n (y) - e 2 ^G (n, s 



d 2 (j) m (y) 
dy 2 



Hn (y) 



2 d A (_ L) r (T) y $n ^n(L)-^ n (-L) 



(18) 



The requirement that the volume variation in (|l8|) vanish yields the SG equations (g), as expected. Of special interest 
to us are surface variations, i.e., the last terms in (|l7|), (|l§|): The requirement that these variations vanish determines 
boundary conditions to (|l|). (For a very clear discussion of the relationship between the surface variation and boundary 
conditions, see Ref. [27| section 11.15.) In the derivation of the surface variation in ([l8|), we have used conditions (^), 
(^) that the local field be continuous at the boundaries y = ±L. The requirement of the continuity of the local field 
at the boundaries x — 0, x = N — 1 has not been so far employed. [Recall our remark in section II that the disregard 
of this requirement is the reason for unphysical solutions to (0), Jl0|).]_ 

If we simply enforce the conditions (|To|), the surface variations in (|17j) and (|lj) cancel out: Thus, all solutions to 
(3), (|ic| ) are stationary points of the Gibbs free-energy functional (|12|). However, under (1C), the surface variation 
in (|18|) does not vanish. Therefore, not all solutions to (^), (|l0|) are stationary points of (13). We have to examine 
conditions of the stationarity of ( jl3| ) in more detail. 

The requirement that the local field is fixed at the boundaries x = 0, x = N — lis equivalent to the requirement 
that the vector potential A is fixed at x = 0, x — N ~ 1. Consider now the total flux In the gauge (fil), we have 



dy [A(N-l,y)-A(fi,y)] 



(19) 



On the other hand, 



$= E J dyH n (y)=HW(N-l) 



N-l 
n=l 



Z7T 



(20) 
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where the first term is the flux penetrating through the boundaries x = 0, x = N — 1, and the second term is the flux 
penetrating through the boundaries y = ±L. (We will call it the "Josephson flux", $j.) Given that 



8A(p,y)=6A(N -l,y) = 0, 



we have 



(21) 



n=l J T y 71=1 



2vr 



= 0. 



(22) 



Thus, the continuity of the field at £ = 0, x — N ~ 1 imposes a constraint on the variations: 

L 



iV— 1 



#n (y) 

dy 



const > 0. 



(23) 



The result ( P3| ) is exactly what had to be expected: By virtue of the Meissner effect, the flux <£>j (and, of course, 
$) is stable against any small perturbations, represented by variations 4> n (y) —* 4> n (y) + 6<j> n (y). Equivalent forms 
of H), H) are 



(L) = 5<f> n (-L) . 



(24) 



4> n (L) - 4> n (-L) = c„ = const > 0. 



(25) 



Note that the existence of conse WW?i|H|^yiii^ a l quantities of the type of $j is a precursor to the existence of soliton 
solutions in nonlinear field theoriesC'lMfl&BEa In Appendix B, we establish a relationship between the conservation of 
$j and the conservation of the intralayer current, which, in turn, establishes a relationship to the variational principle 
for infinite (N = oo) layered superconductorsBlcl 

What will be shown now is that all the stationary points of ([l3]) are the unique minimizers of both (H) and 
First, we notice that high- order variations of ( fl2"| ) and ( |l3| ) coincide: S n F — 5 n £l, n > 2. Thus, all the minimizers of 
( |l3| ) are minimizers of (|12]). On the other hand, the minimizers of (12) obeying (|23| ) are minimizers of (13): From the 
condition of the minimum (Eq), we get 



F 



d(f> n dS(j) n 
dy dy 



H 



F 



dy 



= F 



dy dy 



F 



dy 



H 



> 0. 



(26) 



[Physically, this fact means the equivalence of the description in terms of the Gibbs free energy and the Hclmholtz 



free energy: Because of (E3J) , the functional F 



F 



dy 



;0 



can be regarded as the Helmholtz free- 



energy functional.] Using the standard technique,c3 it is straightforward to prove that all stationary points of (|T^), 
obeying (|23|), are minimizers of jl^ ) [and, thus, of (|13j)]. Indeed, let {<fi n } be the stationary point for $j = $j and 
corresponding H = H. In the vicinity of {</>n}, i-e., for <j) n = <j) n + 6<j> n , we have the following estimate: 



dy 



>r(T) 



N-l N-l 



^ ^G(n,m) / 

1 ™ i «* 



n— 1 m— 1 



d(/> n (y) d(/) m (y) 
dy dy 



N-l 



1 



m V 

n=l 



(/(/ 



d(f> n (y) 
dy 



> -AHr (T) 



(27) 
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Inequalities (27) show that SI has a lower bound in the vicinity of any stationary point {4> n }, obeying (|2^); hence, 
{4> n } is a minimizer of SI and F , Q.E.D. 

To strengthen (p7|) , we minimize the right-hand side of the first inequality with respect to , obtaining 



N-l 



e 2 ^G(n,m) 



m= 1 



dy 



= $n, 

7T 



(28) 



dy 



H 



> -2Hr (T) 



Taking into account that {<^ n } is a solution of (g), making use of ( |ll| ) and (|2q), we get: 



S> 



iff 



> n 



dy ' 



(29) 



r(T) 



N-l N-1N-1 y 

^ — i — i ™ — i « 



#n (y) rf^m (y) 

dy 



iV-l 



> r (T) WJ2[1- cos 0„ (L)] > 0. 



71=1 



(30) 



The inequality SI (^f> n , ; 77^ > is a manifestation of the Meissner effect and had to be expected from general 
thermodynamic arguments. @ Relations (|30|) immediately yield 



dy 



> F 



dy 



N-l 



>r{T)wY^ [1 - cos (i)] + 4ffr (T) $j > 4idr (T) $j. (31) 

3ll) are typical of soliton physics: They are used to establish 



n=l 



We want to emphasize that inequalities of the 
the existence and stability of soliton solutionsJ-M 

Note an alternative interpretation of the variational principle 



sn 



dy 



SF 



dy 



N-l 

n=l J L 



dy 



#n (y) 
dy 



(32) 



The field H in ( |32| ) can be considered as a Lagrange multiplier, implying that variation can be performed without any 
restrictions on 6(j> n (L), 5cj> n (-L). In this case, the requirement that the surface variation vanish yields conditions 
(||); (§)■ Boundary conditions on (p n are uniquely determined by (|23|): see the next section. The value H = H for a 
concrete minimizer $j} should be found from the condition of thermodynamic equilibrium 



an 



H 



0$; 



(33) 



Indeed, SI 



(X 



H ] can be written as 



SI 



dy 



ff =r(T) 



AT-1 

71=1 



1 - cos 0„ (£) - — / dy(j)n (y) sin (f> n (y) 
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4> 



(34) 



which by virtue of (|33|) immediately yields (|9j). 

In summary, we have proved the following: The SG equations (0) and relations (|Io| ) ensure only the stationarity of 
the Gibbs free-energy functional ([jj]) . The necessary and sufficient conditions of the minimum of bot h (O ) and ( |l3| ) 
(which is the Hclmholtz free-energy functional for H = 0) are given by (|j), ||) and the constraint (|23|). Solutions 
to (||), ([l0|) that do not obey this constraint are absolutely unstable. The character of this instability can be easily 
established. Indeed, such solutions are not even stationary points of the Helmholtz free-energy functional, therefore 
S 2 fl — 5 2 F need not have definite sign. Moreover, the functional ([l2]) is unbounded in the vicinity of these solutions. 
Thus, they are nothing but saddle points of (^2|). 



IV. THE PROOF OF THE STABILITY OF THE MEISSNER SOLUTION AND VORTEX-PLANE 

SOLITONS 



A. Boundary conditions on 



Our aim now is to establish boundary conditions on <p n directly from the constraint (p3|). Given that the SG 
equations (|^), boundary conditions (||), (^) and the constraint (|2^) represent a full set of necessary and sufficient 
conditions of the minimum of the Gibbs and Helmholtz free-energy functionals, we will obtain, in this manner, the 
sought proof of the stability of the Meissner solution and soliton vortex-plane solutions. 

First, we observe that since a minimizer of (|l2|), ( |l5| ) must ensure the vanishing of both the surface and volume 
variations in ( |l7| ) and (|l8|), it should necessarily belong to the class of functions that satisfy (||), (||), (§) and the 
symmetry relations (|^). Thus, the variation of the surface terms in ( |l7j ) and ( |l8| ) is performed with respect to trial 
functions that take only discrete values at y = 0: 



4>n (0) = TrZ n , 



(35) 



where Z n can be arbitrarily chosen from the set 0, ±1, ±2, . . . .These functions can be subdivided into classes param- 
eterized by an (N — l)-dimensional "vector" 



Q — (Z\, Z 2 , . . . , Zn-i) ■ 
In view of (j35|), the requirement of the continuity of variations can be met if and only if 

64> n (0) = 0, 



(36) 



(37) 



which means that all the minima of (]l^), ([l3|) are parameterized by the vector Q, and the variation of the surface 
terms in ( |l7j ) and ( |l8|) is performed with respect to trial functions that belong to a certain class (|36|). Moreover, the 
symmetry relations (H) imply 6<fi n (L) = —S(f> n (-L). Combined with (Rh, this yields 



Now we combine dlq) with (p) to obtain 



Scf> n (L) = 8cj> n (-L) = 0. 



4> n (L) - nZ n = — > 0, 



(38) 



(39) 



- 4> n {-L) + ttZ„ = y > 0. 

Since the inequalities in (|39|), ( p0| ) should hold for any Z n , including Z n — 0, we get 

<t>n (-£) < 0, <t>„ (i) > 0. 
Moreover, since for any fixed set {c„} the set {<f> n (±L)} must belong to a certain unique class (^ 

- 7T < 4> n (-L) < 0, 



(40) 



(41) 
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2nZ n < c„ < 2n (Z n + 1) , Z n = 



where [u] is the integer part of u. 

Given (j37|), (38), the boundary conditions ( [5| ) and 



1 f , d(/) n (y) 

-L 



= 0,1,2,. 



(42) 



4> n (0)=7rZ n , Z n = 0,1,2,. 



(43) 



together with (||), (||) determine, in principle, a complete set of conditions for the minimizer of (|l2|), (|l3|). The solution 
belonging to a certain class (|3g) [with Z„ as in ([l3|)l first appears when all </>„ (— Z) = 0, and c n = 27rZ„ [see the 
estimates ( |30| ) , ( |3l| ) : under these condition^-the. density of the Josephson energy at the boundary is a minimum] . 
Thus, we are confronted with the standardEol3l§0 soliton boundary value problem 



n (-L) = O, <M0) = ^Z„, Z„ = 0,l,2, 



(44) 



plus the .boundary conditions (^) . To find out what type of minimizers can be realized in reality, we have to solve 
Eqs. (^).lO Note that N — 2 relations (|]) implicitly impose N — 2 conditions on N — 1 constants Z n . Therefore, as 
shown in Refs. ||,[l7|, the only solutions to (||), J44|), compatible with (||), are those that satisfy (Q), i.e., the Meissner 
solution and the soliton vortex-plane solutions: 



= (N V ,N V ,...,N V ). 



1 

27 



dy 



dy 



0,1,2,. 



(45) 



The progecties of these solutions will be discussed in more detail in what follows. Here we want to brief on the results 
for H = 0.EO For H = 0, L < oo, there are no soliton solutions at all (including the vortex planes), and the only stable 
solution is the trivial Meissner solution <pi = <p2 = ■ ■ ■ 4>n-i = 0. The situation changes drastically for H — 0, L = oo: 
The imposition of soliton boundary conditions on (j) n (±oo) automatically ensures the fulfillment of the boundary 
conditions ^g 5 - (±oo) = by virtue of Eqs. (^) themselves and some elementary theorems of mathematical analysis. 
Aside from the vortex-plane solution with N v = 1 in ( |45| ) (where L = oo), we have a variety of soliton solutions ( |3^ ) 
with Z n arbitrarily chosen from the set 0, ±1. The fact that for -ffrf= 0, L — oo each <j) n can "accommodate" no more 
than one vortex or anivortex is a generalization of the well-known&El result for the single static SG equation and can 
be easily proved by the use of the first integral (O) with the right-hand side equal to N — 1. 



B. The Meissner solution and soliton vortex-plane solutions 



The range of the existence of the solutions, parameterized by (f45|), is determined by the boundary value problem 
(|l|), (H) (with all Z n = N v ) and relation ©lO 



0<H<H = H sL , for N v = 0, 



(46) 



m s <H<H N ^ for N v = 1,2, 



(47) 



where H s l has the meaning of the superheating field of the Meissner state (N v — 0) for L < oo (H s l > H s for 
L < oo, and H soo = H s ). The lower bound in (Jl6|), p7| ) is determined by the exact upper bound in (|4l]). At 
H = supH = Hn v , when all 4>n = m f 4>n {~L) = ~ tt> there is instability of the saddle-point type (see the end 

of section III). Note that both the Meissner solution [N v = 0) and the vortex-plane solutions (N v > 1) automatically 
satisfy the symmetry relations (j^).SO 

It is instructive to verify the general inequalities (^o|), (|3l|). Mathematically, it is sufficient to do this only for H 
equal to the lower bounds in (pl6|), ({47]): By continuity arguments, the result will be valid in the whole field range. 
For the Meissner solution the verification is trivial. For N v > 1, we employ the exact expression 



F 



dy 



H 



N-l N-l 



( T ) g2 E E G ("' m ) / d y 



n—l m— 1 



d4>„ (y) dcjhn jy) 
dy dy 



(48) 
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where H = , / Hfj _ 1 — if| . As shown in Ref. in this case 



d<t> n (y) 

dy 



= ^(±L) = 2y/H% v _ 1 -H*, for all n. (49) 



Combining Eqs. fl48|), (j49), we get exactly the lower bound in (|3l|). 

As is emphasized in Ref. |[ the obtained solutions are valid for any N, including the cases N = 2 (a single 
junction) and N = 3 (a double-junction stack). (For N — 2, 3, we have derived in Ref. |^ exact, closed-form analytical 

expressions.) The solutions with N v > 1 are pure solitons only at H — \J Hfj ~ x — H^, when j n (±£) = sin^> n (±£) = 

[j n (y) is the density of the Josephson current© for /„ = 1] . In the rest of the regions_(E7]) , we have solitons " dressed" 
by the Meissner field. In the case of N — 2 (the single junction), Owen and ScalapinoEScalled these regions the "N v to 
N v + 1 vortex mode" . [Because the principle of superposition does not apply to the nonlinear Eqs. (|2|), the Meissner 
and the vortex-plane fields cannot be separated from each other.] It is clear that the vortex-plane solutions for N > 2 
are a direct generalization of ordinary vortices in single junctions. 



Of special interest is the overlap of the regions (M), (47) for N v = N v and N v = N v + 1. As a result, the obtained 



solutions cover the whole field range < H < oo, as they should. Mathematically, the overlap is related to the fact 
that the solution with N v = N v cannot be continuously transformed into the solution with N v = N v + 1 by changing 
H, as is always the case for solitons. For the single junction, the overlap was first established numerically in Ref. |3j 
and discussed qualitatively in Ref. The overlap practically vanishes for Hff v ~> H s . Given that all H^ v decrease 
when W = 2L increasesjj the overlap is stronger for large W and can involve several neighboring states. Physically, 
the actual equilibrium state is the one that corresponds to the absolute minimum of the Gibbs free energy for given 
H. The rest of the allowed states are metastable. In view of the above-mentioned discontinuity, a transition from 
the state with N v = A^-to-the state with lower Gibbs free energy N v = N v + 1 will necessarily be a phase transition 
of the first-order type.EJoEa It particular, the lower critical field H c \ is determined from the requirement that the 
Gibbs free energy of the state N v — 1 be equal to that of the Meissner state (N v = 0) and satisfies the relation 

^Hj L - m s < H cl < H sL . 



C. Topological considerations and stability in the dynamic regime 

The stability of the Meissner solution and vortex-plane solitons can be better understood, if we analy^-tie-nbtained 
results from the general point of view of the stability of topological defects in continuum media. □otffitjSEj To this 
end, we consider the density of the Gibbs free energy ( |l2| ) at the boundaries y_ = ±L. 

Because of the general symmetry relations (g), valid for any solution to (|6|), (|o|), the density of the Josephson 
energy is equal at y = —L and y = +L: 

1 — COS (j> n (— L) = 1 — COS (f) n ( + L) , n — 1 < X < U, Tl = 1, 2, . . . , N — 1. (50) 



Taking into account (|l l|) , we conclude that also the density of the total free energy is equal at the boundaries y = —L 
and y — +L and thus corresponds to the degenerate equilibrium ("vacuum") state, unperturbed by topological defects 
(solitons). Mathematically, the boundary of the interval — L < y < L can be considered as a 0-dimensional sphere: 
S° = {—L, +L}. Given that configurations <fi n and cj) n + 2irZ n (Z n = 0, ±1, ±2, . . .) are physically indistinguishable, 
we can fix the values (— L) as in (^d]) and regard the functions 

cj) n (+L) + dp n {-L) 

ip n (+L) = = Z n (51) 

as continuous maps of the boundary into the additive group of the integers, Z: S° ^> Z. (Z is the group of the 
degeneracy of the equilibrium state, or the order-parameter space.) The fact of the existence of topological! 
maps of this type, realized by soliton solutions, is often expressed in terms of the "zeroth homotopy group"! 
ttq (M), where the index "0" stands for the boundary 5° and M is the order-parameter space: 

tt (Z) = Z. (52) 

[Note that ttq (M) is merely the set of disconnected components of the space M .] Because of the boundary conditions 
([|), all tp n at H > realize the same mapping: Zi = Z 2 = ... = Z^-i = Z. The external field H > breaks the 
symmetry <f> n — > — <p n [see the second term in (U3)]. Therefore, only the values Z = N v = 0, 1, 2, . . . are allowed, with 



N v — being the "vacuum", Meissner state. In this way, we arrive at the natural topological classification (45) of the 



11 



minimizers of ([l2]). Owing to the continuity conditions (|36"|), (|37|), variation in (p7|), ( |lli| ) is allowed only with respect 
to trial functions {4> n } that have the same end points <f> n (±L) and the middle point <fi n (0) as the minimizer |0 n }, 
i.e., {4> n } are homotopic to {4> n } and belong to the same class ([45]): hence the stability of {(f> n } against continuous 
perturbations. 

Numerical simulations for time-dependent coupled SG equations have revealed exceptional stability of vortex planes 
in the dy namic regime as well: see figure 7 in Ref. [50| figure 8 in Ref. |3l], and figure 3 in Ref. |3^. (The authors of 
Refs. |30[f3l|]3"^ employ the terms "coherent", "in-phase" or "phase-locked modes" instead of our term "vortex planes" 
that we prefer for physical reasons.) Although a detailed analysis of the dynamics of vortex planes is beyond the 
scope of this paper and will be done elsewhere, the results of Refs. 30 , 3l|,^^ can be explained already at this stage. 



In the absence of dissipation, the dynamic SG equations, describing an evolution of the system in the time interval 
ti < t < tf, can be derived from a corresponding Lagrangian by use of a variational principle. The requirement that 
the surface variation vanish on the whole perimeter of the space-time boundary leads to the conditions 8<j) n {y,ti) = 
5(j) n (y,tf) — and a generalization of the conservation law for the Josephson flux <£>j, Eq. (p3|): 



N-i h ^ , / ,x N-l 



rh V" a 1 f a d( t )n *) rK <t>n(L,t) -<j) n (-L,t) . , 

= 2^ *«2^ / d y — qt~ = 1^ $ « ^ = cons 1 (53) 

n=l _ L y n=l 



which means that the differences <j> n (L, t) — cf> n {—L, t) do not depend on t. Thus, by fixing the boundary conditions 
4> n (±L, t) at t = U as in subsection IV. A, we fix the initial value of the flux $j = $j that will not change in the 
course of the evolution of the system from t = ti to t = tf. The topological type of the solution [see (^)] will not 
change, either: 



N v = 



J_ f dy d^(y,t) 



2n J dy 



-L 



4>n (L,t) - (f) n (-L,t) 

2tt 



= const. (54) 



As usualJjQ'LlEHj this situation can be formalized in terms of the conserved topological current 

AT-l « 

3ii = 2_ \ TT^^d-iK; df.jf, = 0, (55) 

n=l 

where fj,,j = 0,1; 9 M = (dt,d x ); and e MT is the antisymmetric symbol on two indices, eoi = — eio = lj with $,/ = 

L 

f dyj being the topological charge. 

-L 

As should be clear from these results, time-dependent SG equations alone cannot describe Josephson-vortex pene- 
tration, i.e., an evolution of the system from the topologically trivial Meissner state, N v — 0, to a state with N v ^ 0. 
Unfortunately, this important issue has not been realized in Ref. [l2] that claims to have " demonstrated a dynamic 
process of vortex penetration" by means of numerical simulations for time-dependent SG equations. 

V. DISCUSSION 

A. Unstable solutions to the SG equations 

As is proved in section III, all non-topological, non-soliton solutions to (@), ([l(]) that do not meet the requirement 
of the flux conservation (^3|) are absolutely unstable: They are nothing but saddle points of the Gibbs free-energy 
functional (|l2|), cannot be assigned any "free energy" and are therefore unobservable. Since the requirement of the 
continuity of variations in ( |l7| ) docs not impose on such solutions any constraints of the type (36), (pTj), they can 



be continuously transformed into the stable Meissner solution or a vortex-plane solution, representing the actual 
minimum o f (|l2|) at a given H, by a series of infinitesimal deformations of <p n without a violation of the boundary 
conditions (|10[) . [A clear illustration of such a transformation for non-topological defects in a system of planar spins 
see in Ref. [29f section II. B.] 

Analytically, all unstable configurations for H > can be obtained using the symmetry relations (^), as solutions 
to the boundary value problem 

^(-L) = 2H, (b n {0) = irZ n , Z n = 0,1,2,... (56) 
dy 
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that violate topological boundary conditions, derived in subsection IV. A. For example, one can set in ( pq ) Z\ = Z2 = 
. . . = Z N _x = Z and increase H beyond the upper bound H Nv of the stability regions jj^), Jl7| ) for a given = 
By continuously increasing H beyond the stability region of the Meissner state, unstable configurations with Z = 0, 
interpreted as " Josephson- vortex penetration", were obtained in numerical simulations for static SG equations (Ref. 
|3l| , figure 7) and time-dependent SG equations (Ref. O, figure 2). Analogous instability for Z > 1 is demonstrated 
by numerical simulations in the dynamic regime [H =const > 0, and an increasing transport current I > I c (H)] in 
Ref. |3^: see region IV in figure 3 therein; region III corresponds to dynamically stable vortex planes. 

Unstable solutions appear also when not all Z n in (|5^ ) are equal to each other. Thus, an unstable "single Josephson 
vortex" corp-sjinnds to the choice Z\ = 1, Z n= ^\ — 0. Solutions of this type were obtained in several numerical 
simulations .ESoEj By way of illustration, wc consider here only the case H = 0. As is explained at the end of 
subsection IV. A, the only stable configuration for H — 0, L < 00 is the trivial Meissner state 4>x = <j>2 = ■ ■ ■ = 
4>n-i = 0. Figure 5 in Ref. and figures 1, 2 in Ref. [if] clearly show that the solutions presented therein, in 



reality, are characterized by all Z n = 



-L 



d<t> n (y) 
dy 



= and, thus, belong to the class Q„ = (0, 0, . . . , 0) of the 



general topological classification (|45|). By means of continuous deformations, they can be transformed into the trivial 
Meissner solution. 

Other unstable solutions, available in the literature, can be analyzed along the same lines. In particular, the 
"triangular Josephson- vortex lattice with the period x p = 1", proposed in Ref. [19, corresponds to the case Z dd = Z , 

^even Z 



B. A comparison with Abrikosov vortices in type-II superconductors 

As the formation of a vortex-plane soliton involves only phase differences between successive S-layers, it does 
not affect the topology of the layered superconductor. In contrast, the appearance of a linear (R 1 ) singularity of 
the order parameter A (r) = |A (r) exp [iip (r)] is necessary for the formation of an Abrikosov vortex in continuum 
type-II superconductors. Thus, in the presence of a single Abrikosov vortex, the topology of the continuum type-II 
superconductor changes from R 3 (the three-dimensional Euclidean space) to R3/R 1 = S 1 . Therefore, the notion of 
the "vortex core"llfl is inherent (both physically and mathematically) to the Abrikosov vortex and is meaningless in 
the case of the vortex plane. 

An isolated Abrikosov vortex is itself a stable object, both topologically and energetically. (The latter can be proved 
by the use of the same mathematical methods as those employed in our section III: see, e.g., Refs. (HI) Therefore, an 
equilibrium state of N v Abrikosov vortices is determined by comparing the values of the Ginzburg-Landau free-energy 
functional for different spatial configurations, which yields the well-known triangular lattice as the most favorable 
one.EJ In contrast, the notion of the " Josephson- vortex lattice" is senseless for layered superconductors: One can only 
speak of A t ,-soliton (vortex-plane) states, with N v = representing the Meissner state, and each vortex plane being 
a " Josephson vortex" itself. 

In the case of extreme type-II superconductors, the linear singularities, associated with Abrikosov vortices, can be 
easily incorporated into the framework of the simple London model.Lj The resulting equation is a linear inhomogeneous 
partial differential equation for the local field. Owing to linearity, the local field is a superposition of the Meissner 
and vortex fields. As is emphasized in subsection IV. B, this is not the case for layered superconductors because of the 
nonlinearity of the SG equations (|^). Unfortunately, this important issue was not understood in some publications 
concerned with Josephson-vortex penetration.! 3 ^ 

As we can see, there is no " analogy" between the Abrikosov- vortex structure in continuum type-II superconductors 
and the Josephson-vortex structure in layered superconductors in the naive sense. til Instead, there is a much subtler 
mathematical analogy: The topological classification of vortex configurations in type-II superconductors is isomorphic 
to that in layered superconductors. A proof is straightforward. For the reasons explained above, the boundary of a 
type-II superconductor is, in general, topologically equivalent to a one-dimensional sphere (a circle) S . The order 
parameter space is M — U (1) (the symmetry group of quantum electrodynamics). Topologically, U (1) = S . Thus, 
soliton solutions, in this case, realize nontrivial maps S 1 — ► S , l in A4J r] the-|Continuous maps S 1 — > S 1 have a group 
structure of the fundamental (or first homotopy) group tti (S* 1 ) £I'EIeH3qEIE3 Given that S 1 = R/Z (R is the additive 
group of the real numbers), we can write 

7T! (S 1 ) = 7T! (R/Z) = n Q (Z) = Z, (57) 

which should be compared with (|5^). As in the case of layered superconductors, the external magnetic field H > 
breaks the symmetry tp — > —ip. Thus, only the states parameterized by N v = 0, 1, . . . are possible, with N v = being 
the "vacuum", Meissner state. 
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To conclude this discussion, we have to clarify a typical misunderstandingeHl concerning the role of the laminar 
modelEJ in type-II superconductivity. In reality, the order-parameter space of a continuum type-II superconductor, 
M = S x , precludes the existence of topologically stable plane defects, envisaged by the laminar model. Indeed, 
consider two points P\ = (xo,a, zq), Pi = (xo,b, zq) on the opposite sides of such a defect, in unperturbed regions 
of the superconductor. Join these points by a continuous path, parameterized by a < y < b. The boundary of the 
interval [a, b] is a O-dimensional sphere S° — {a, Z^juhich. leads us to a consideration of the maps S° — > S . However, 
the pertinent homotopy group ttq (S 11 ) is trivial ,Btffli§Bc3 i.e., 

{S 1 ) = 0, (58) 

in contrast to ( |5^ ) and (|57|). Topological instability of the "laminar solution" hardly allows one to expect that this 
solution corresponds to any local minimum of the Ginzburg-Landau free-energy functional. Therefore, a comparison 
with unstable "isolated fluxons",-er "triangular Josephson- vortex lattices" in layered superconductors is much more 
appropriate than the far-fetchedt3 "similarity" to the vortex plane. 



C. The interpretation of experimental data 

Experimental observations of the vortex structure in layered superconductors-tpjiJpe roughly subdivided into two 
groups: (i) direct observations, allowing oneJtnllida.ualize" the fluXpdistributionjE^lnJO (ii) indirect observations (i.e., 
measurements of c-axis transport properties,tj'LJL3 magnetization,o and the upper critical fieldcil). Here, we present 
an overview of these observations, showing that all the experimental data available up to now can be explained in 
terms of the stable vortex-plane configurations. A detailed quantitative analysis can be done with the use of the 
results of Refs. g,|§ 

Josephson-flux distribution, characteristic of vortex planes for H > 0, was directly observed on artificial low-T c 
stacked junctions in Ref. ^ (by low-temperature scaaning electron microscopy) and in Ref. |^| (by polarized neutron 
reflection). In particular, the double-junction stacked (N = 3) has revealed the phase-difference symmetry (fix = fo, 
exactly as could be expected from the general relations (||) for vortex planes. In both the experiments,E3o penetration 
of the flux occurred simultaneously and coherently into all the junctions, in full a gree ment with the scenario for the 
vortex planes.Era Moreover, accompanying measurements of magnetization in Ref. |24| have shown typical oscillations 
and hysteresis. The oscillations should be viewed as a manifestation of a series of first-order phase transitions, 



discussed in section IV. 2 and Refs. [l||2|, whereas the hysteresis is implied by the overlap of the regions fl46|), (|47|). 

We draw attention to a possible application to high-T c superconductivity. Oscillations of magnetization in parallel 
fields, interpreted as evidence of Josephson nature of the flux, have been reported for YBCO in Ref. According 
to Ref. "the temperature dependence of thermagnetization contradicts the present theoretical expectations". 

As shown by our self-consistent calculations,™ 2 ! the oscillating behavior of the critical Josephson current I c (H) 
(the Fraunhofer pattern) is a result of successive penetration of vortex planes and their pinning by the edges of— the. 
superconductor. Oscillating I c (H) dependencies-have been observed both on artificial low-31-^tacked junctions' 2 ^ " v 



and high-T c layered superconductors BSCC0.E3O "Irregularities" of the dependence I c (H)f3x3 such as, e.g., mul 



ap of the regions (phf), (f47|). Behavior of this type 
which confirmed the theoretical prediction of the 



tivaluedness and aperiodicity, can be easily explained by the ove 
was observed a long time ago on the single Josephson junction,^ 
overlap for ordinary Josephson vortices.tfl 

The most " ancient" experimental confirmation of the stability of jrortex planes is provided by observations of the 
"crossover" behavior of H C 2 (T) in artificial low-T c stacked junctions.cil For H < H C 2, the condition f n = 1, employed 
in the derivation of Eqs. (||), in no longer valid. However, periodic modulations of /„ (y), caused by the presence of 
vortex planes and therefore identical in all the S-layers, account for the observed behavior of H C 2 (T).Era 

Finally, we want to-comment on direct observations of non-equilibrium isolated vortices in layered high-T c super- 
conductors at H = 0H3 As is explained at the end of subsection IV. A, nontrivial flux configurations cannot exist in 
a layered superconductor with ideal periodicity at H = 0, L < oo. However, the presence of structural defects (e.g., 
stacking faults, as is hinted in Ref. [37]) violates the condition of ideal periodicity and should stabilize energetically 
an otherwise unstable configuration. In this situation, we indeed expect to obtain non-equilibrium isolated vortices, 
because their self-energy is lower and the c-axis extent is smaller than those of vortex planes at H — 0, L < oo.El A 
detailed mathematical analysis of this case can be done on the basis of the results of Ref. 0. 



VI. CONCLUSIONS 

In brief, we have solved the problem of the classification of all solutions to (||), ([l0]) with respect to their stability. In 
our consideration, we have employed exact methods of the calculus of variations, soliton physics and the exact results 
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of Ref . H [the expression for the Gibbs free-energy functional ( [l2| ) , the first integral (|Tl|) and the solution to the soliton 
boundary value problem (|J), (Q), (0)]. In view of obvious mathematical complexity, the problem of the stability of 
vortex configurations could not be solved by use of inadequate methods, employed, e.g., in Refs. |ll 12 T^, [l9| , p0pl| , 
and was not even posed in any of these publications. 

In full agreement with the fundamentals of soliton physics, we have proved that the only minimizers of both the 
Gibbs and Hclmholtz free energy functionals are the Meissner solution (the "vacuum" state) and soliton vortex-plane 
solutions. They represent the actual equilibrium field configurations. The obtained results allowed us to explain 
exceptional stability of vortex planes, established in numerical simulations, and to provide a unified interpretation of 
the experimental data available up to now. 

In contrast, non-soliton configurations ("isolated fluxons", triangular Josephson-vortex lattices", etc.), proposed in 
previous publications, turned out to be absolutely unstable and unobservable: They are nothing but saddle points of 
the Gibbs free-energy functional and are not even stationary points of the Helmholtz free-energy functional. Physically, 
these configurations violate conservation laws for the flux and the current. 

One may ask a natural question why exactly the unstable configurations were previously proposed as the " equilib- 
rium state" , whereas the actual minimizers of the free energy (vortex-plane soiitons) were neglected. We think that 
the answer lies in the following: 

i) the hypothesis of an "analogy" with the Abrikosov- vortex structure in type-II superconductors, accepted without 
any mathematical justification^ 

ii) the absence of an exact mathematical definition of the "Josephson vortex". In Refs. [ll 18, "isolated Josephson 



vortices" were discussed without any consideration of the SG equations. The fact that Josephson vortices are nothing 
but static soiitons of the SG equations was not realized in subsequent publications, either. For example, the existence 
of Josephson vortices in the case W <C 2Aj was denied in Refs. [l5|,|l6|. which is refuted by our Eqs. (|47|), valid for any 

w :.. n 

iii) obvious mathematical mistake s in the treatment of the Lawrence-Doniach modell3 for infinite (N = oo) layered 
superconductors. The neglectt-3il3'E3Ell of the surface variation in the variational principle for the Lawrence-Doniach 
functional resulted in a loss of the conservation laws for the current and the flux, as was first pointed out in Ref. EE: 
Soliton solutions are a corollary of these conservation laws; 

iv) the absence of any investigation of analytical properties of the coupled static SG equations for H > 0, W < oo. 
Pertinent soliton solutions were obtained in our papers: in the exactly solvable cases N = oo ffirfrjJKp N — 2,3 
(Refs. 10), and in the general case 2 < N < oo (Ref. |J). Standard methods of soliton physicsDfmQmJH as well as 
advanced methods of the calculus of variations and of the theory o£ HifWpntial equations, employed in our analysis, 
were completely disregarded in previous theoretical publications Jl3'Ej'E3til 

v) the alispnrp of any ^t.t.pmpt.s to analyze the stabili,^ of the proposed "vortex" configurations, both in theoret- 



ical publications 



ii 
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and numerical simulations .113 Dynamic stability of vortex planes, noticed in numerical 
simulations ,E3EilS was not understood and neglected-Hi— ■ 

vi) the neglect of direct experimental observationscSEJ of the Josephson-vortex structure at H > 0: These obser- 
vations have clearly revealed that exactly the vortex planes (not "isolated fluxons" or "triangular lattices") are the 
actual equilibrium field configurations; _ . 

vii) long-term doiaaination of the subjective point of viewEU and the absence of any pluralism of opinion. As a result, 
the critical remarksca am-neglected, whereas the attempts to clarify the situation within the framework of a rigorous 
mathematical approaclffil are immediately attackec£3 with the use of inappropriate methods £3 

We hope that this paper will finally convince both theorists and experimentalists, specializing in the field of 
weak superconductivity, of the necessity to give up the old, unsound theoretical prejudices: The wealth of magnetic 
properties of layered superconductors (both low- and high-T c ) cannot be understood without the soiitons. One should 
also think of possible practical applications of the vortex-plane soiitons, e.g., in submillimeter-wave generators,.-as is 
proposed in Ref. |32|. Given the role of the single SG equation indifferent fields of physics (quantum optics;Q the 
Skyrme and the Thirring models in elementary particle^hysicsjtffiO the theory of dislocations and magnetism, let 
alone the Josephson effect, in condensed matter physics,Q etc.), we expect that our exact results for the coupled SG 
equations (including the single one as a particular case) may find applications in these fields as well. 

The coupled SG equations for H > 0, W < oo have not been studied in mathematical literature, either. Our exact 
analytical results for the static case constitute only the first step in this direction. The next stage should be analytical 
properties of time- dependent equations. Our paper may stimulate interest in this problem of specialists in applied 
mathematics as well. 
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APPENDIX A: THE PROPERTIES OF THE MATRIX G (n, m) 

The explicit form of G(n, m) is 

, ,n ( ..m—N ,.N—Tn\ _i_ ,.N — n — m ,.m\ 
\n-m\ _ f 1 {f 1 ) + V Ip ~ P 7 



G(n, to) 



2ejl 



where 



e 2 / ? 
p = 1 + T -'V 1+ 4- 



The following properties of G(n,m) are obvious: 

G{n,m) = G(m,n), 

G(n, N — m) = G(N - n, m). 
The matrix G(n, m) is positive definite, since all its eigenvalues are positive: 



A 2 

e,- = -4, A., = 



-J £ 2 



2 + e 2 -2cosff 



Of importance are the summation rules: 

JV-l 



./ 1-2 V 1. 



G(n, to) = ^ t 1 - G{n, 1) - G{n, N - 1)] 



(Al) 

(A2) 

(A3) 
(A4) 

(A5) 



m— 1 



l - 



1< n < N- 1: 



(A6) 



where 



JV-l JV-l 



1] H G(n,TO) = i 



n— 1 m— 1 



N - 1 - 



2Jl + ^- ei _^-i 



TV — 1 



V^T- 6 !-^- 1 
e(iV-l) l + n N 



(A7) 



(N — 1)N 



f— l-i 

.2 _ 2 7T(2fc + l) 



fc=0 



is the superheating (penetration) field of a semiinfinite (0 < y < oo) layered superconductor! 



(A8) 



APPENDIX B: A RELATIONSHIP TO THE VARIATIONAL PRINCIPLE FOR INFINITE LAYERED 

SUPERCONDUCTORS (N = oo) 

The intralayer currents for /„ = 1 are given b; 

ay 



— [H n (y) - H n+1 (y)} , n = 0, 1, . . . , N - 1, 



where Hq (y) = Hn (y) = H. Using the second relation in (Bl), we get 



N-l 



n=0 



which is the conservation law for the total intralayer current. Moreover, in view of (§ and Q, (|A4|), 

H n (y) = H N - n (y) . 

Hence, 

Jn (y) = -Jn-ti-x (y) ■ 



Using the first relation in (Bl), we write 



dy[J n (y) - J n -i (y)] 
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(Bl) 

(B2) 

(B3) 
(B4) 



4> n (L) - <f> n (-L) - I dy[A(n,y)-A(n-l,y)}, n = 1, 2, . . . , TV - 1. 

-L 



(B5) 



The second term on the right-hand side of (B5) is the flux between the S-layers n and n — 1. We can therefore rewrite 
© using (|): 



dy [J n (y) - J n _x (J/)] = W [G (n, 1) + G (n, JV - 1)] 



JV-l 



6„ (L) - K (-L) + - E G (n, m) [0 m (I) - m (-1)] , n = 1, 2, . . . , TV - 1. 



In view of the flux-conservation conditions (24), the variation of the right-hand side of (pq) vanishes, hence 



6J n (y) = 6J n -x{y), n = 1, 2, . . . ,N - 1. 



Combined with the current-conservation law flB2|), relations (B7) yield: 

«J„(») = 0, n = 0,1,..., JV-l, 



(B6) 



(B7) 



(B8) 



which means that partial intralayer currents are also conserved. 

Consider the case JV 3> 1. For n satisfying the condition [e _1 ] <C n <C N — 1 — [e _1 ] , we can proceed to the limit 
N — > co in the second relation (Bl), obtaining 



J n (y) = ^-2A(n,y)=0. 
dy 



(B9) 



This is exactly the result derived for the infinite (JV = co) layered superconductor in Refs. [l],^| by means of an exact 
variational principle, based on the use of the conservation law for the total intralayer current. 
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